
# set up the working directory in your computer
setwd("")

# ---------------  figure 1 
data <- read.csv("OECD_fertility.csv")
data <- subset(data, Year!=2014)

year <- data$Year
country <- names(data)[-c(1,ncol(data))]

pdf("figure1.pdf",height=5,width=6,paper="special")

par(mar=c(5.1,4.1,2.1,2.1))
plot(year,year, type="n", ylim=c(1,8), ylab="Total Fertility Rate(TFR)", xlab="Year")
for (i in 2:ncol(data)) if (i != ncol(data)) lines(year,data[,i], col="gray", lty=3)
lines(year, data[,which(country %in% "KOR")+1],col="black", lwd=2, lty=1)
lines(year, data[,which(country %in% "JPN")+1],col="black",lwd=2, lty=5)
legend("topright",c("South Korea","Japan","Other countries"), col=c("black","black","gray"), lty=c(1,5,3), lwd=c(2,2,1), bty="n")
dev.off()



# ---------------  figure 2
all.data <- read.csv("fertility_descriptive1.csv")
all.data <- subset(all.data, r_age <= 40)
mean.b <- array(NA,dim=c(6,3))

mean.b[1,1] <- mean(all.data$ever_b[all.data$r_work==0&all.data$strong==1],na.rm=TRUE)
mean.b[2,1] <- mean(all.data$ever_b[all.data$r_work==1&all.data$log_rp< 0&all.data$strong==1],na.rm=TRUE)
mean.b[3,1] <- mean(all.data$ever_b[all.data$r_work==1&all.data$log_rp>=0&all.data$strong==1],na.rm=TRUE)
mean.b[4,1] <- mean(all.data$ever_b[all.data$r_work==0&all.data$strong==0],na.rm=TRUE)
mean.b[5,1] <- mean(all.data$ever_b[all.data$r_work==1&all.data$log_rp< 0&all.data$strong==0],na.rm=TRUE)
mean.b[6,1] <- mean(all.data$ever_b[all.data$r_work==1&all.data$log_rp>=0&all.data$strong==0],na.rm=TRUE)


mean.b[1,2] <- sd(all.data$ever_b[all.data$r_work==0&all.data$strong==1],na.rm=TRUE)
mean.b[2,2] <- sd(all.data$ever_b[all.data$r_work==1&all.data$log_rp< 0&all.data$strong==1],na.rm=TRUE)
mean.b[3,2] <- sd(all.data$ever_b[all.data$r_work==1&all.data$log_rp>=0&all.data$strong==1],na.rm=TRUE)
mean.b[4,2] <- sd(all.data$ever_b[all.data$r_work==0&all.data$strong==0],na.rm=TRUE)
mean.b[5,2] <- sd(all.data$ever_b[all.data$r_work==1&all.data$log_rp< 0&all.data$strong==0],na.rm=TRUE)
mean.b[6,2] <- sd(all.data$ever_b[all.data$r_work==1&all.data$log_rp>=0&all.data$strong==0],na.rm=TRUE)

mean.b[1,3] <- sum(!is.na(all.data$r_work==0&all.data$strong==1))
mean.b[2,3] <- sum(!is.na(all.data$r_work==1&all.data$log_rp< 0&all.data$strong==1))
mean.b[3,3] <- sum(!is.na(all.data$r_work==1&all.data$log_rp>=0&all.data$strong==1))
mean.b[4,3] <- sum(!is.na(all.data$r_work==0&all.data$strong==0))
mean.b[5,3] <- sum(!is.na(all.data$r_work==1&all.data$log_rp< 0&all.data$strong==0))
mean.b[6,3] <- sum(!is.na(all.data$r_work==1&all.data$log_rp>=0&all.data$strong==0))

lci <- mean.b[,1]-mean.b[,2]/sqrt(mean.b[,3])*1.96
uci <- mean.b[,1]+mean.b[,2]/sqrt(mean.b[,3])*1.96

out.b <- array(0,dim=c(3,2))
out.b[1,1] <- mean.b[1,1]
out.b[2,1] <- mean.b[2,1]
out.b[3,1] <- mean.b[3,1]
out.b[1,2] <- mean.b[4,1]
out.b[2,2] <- mean.b[5,1]
out.b[3,2] <- mean.b[6,1]

pdf("figure2.pdf",height=5,width=6, paper="special")
barplot(out.b, beside=TRUE, col=c("white","gray","black"), ylim=c(0,0.25), 
	legend=c("Not Working","Log Wage Ratio < 0","Log Wage Ratio >= 0"), args.legend=c(bty="n"),
	names.arg=c("Strongly Embedded Couples","Weakly Embedded Couples"),
	ylab="Fertility Rate")
dev.off()


# ---------------  figure 3
data <- read.csv("predicted_coef.csv")
x <- read.csv("predicted_at.csv")
x <- x[,1]
n <- length(x)

wdata <- data[rep(c(0,1),n)==0,]
sdata <- data[rep(c(0,1),n)==1,]

w.y <- wdata[,1]
s.y <- sdata[,1]


pdf("figure3.pdf",height=6,width=6, paper="special")
plot(all.data$log_rp_hat, jitter(all.data$ever_b,0.2), type="n", ylab="Predicted Fertility Rate",xlab="Log Wage Ratio [imputed]")
points(all.data$log_rp_hat, jitter(all.data$ever_b,0.3), cex=0.8,pch=19,col=ifelse(all.data$strong==0,"black","gray"))
points(x,s.y, type="l", col="gray", lwd=3, lty=2)
points(x,w.y, type="l", col="black", lwd=3, lty=1)
legend("right",pch=19,col=c("gray","black"), lty=c(2,1), c("Strongly Embedded Couples","Weakly Embedded Couples"),bty="n")
dev.off()

